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We explore the effect of a cubic correction <?nl</> 3 on the mass function and bias of dark matter 
0^ , haloes extracted from a series of large N-body simulations and compare it to theoretical predictions. 

Such cubic terms can be motivated in scenarios like the curvaton model, in which a large cubic 
correction can be produced while simultaneously keeping the quadratic /nl0 2 correction small. The 
deviation from the Gaussian halo mass function is in reasonable agreement with the theoretical 
^ • predictions. The scale-dependent bias correction Ab K (k, <7nl) measured from the auto- and cross- 

power spectrum of haloes, is similar to the correction in /nl models, but the amplitude is lower than 
. theoretical expectations. Using the compilation of LSS data in Slosar et al. [JCAP, 08, 031 (2008)], 

we obtain for the first time a limit on # NL of -3.5 x 10 5 < g NL < +8.2 x 10 s (at 95% CL). This 
' limit will improve with the future LSS data by 1-2 orders of magnitude, which should test many of 

the scenarios of this type. 
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I. INTRODUCTION 

In standard single field inflation, primordial curvature 
perturbations are produced by the inflaton field as it 
slowly rolls down its potential ( 0, H H, 0] ) ■ Most of these 
models predict a nearly scale-invariant spectrum of adi- 
abatic curvature fluctuations in agreement with cosmo- 
logical observations. In addition, very small deviations 
from Gaussianity are expected p, H 0] • Therefore, any 
evidence for or against the detection of primordial non- 
Gaussianity would strongly constrain inflationary scenar- 
ios. 

Non-Gaussianity can be produced by nonlinearities in 
the relation between the primordial curvature perturba- 
tion $ (Here and henceforth, the usual Bardeen potential 
in matter-dominated era) and the inflaton field, interac- 
tions of scalar fields, a modified dispersion relation or a 
departure from the natural adiabatic vacuum state (see 
P for a review). Any non-Gaussianity that is generated 
outside the horizon induces a three-point function (or bis- 
pectrum) B$(ki, k2, ka) that is peaked on squeezed tri- 
angles (i.e. k\ <C &2 ~ ks) for realistic values of the scalar 
spectral index. The resulting non-Gaussianity depends 
only on the local value $(x) of the Bardeen's curvature 
potential and can thus be conveniently parametrised up 
to third order by 

*(x) = 0(x) + /nl [^ 2 (x) - {4> 2 )] + 5 nl0 3 (x) , (1) 

where </>(x) is an isotropic Gaussian random field and 
/nl, <?nl arc dimcnsionless, phenomenological param- 
eters. While the quadratic term generates the irre- 



ducible three-point function or bispectrum at leading 
order, the cubic term does so for the irreducible four- 
point function or trispectrum. These statistics can be 
computed straightforwardly from a perturbative expan- 
sion of the homogeneous Robertson- Walker background 
d [l~fj| . Convolved with the appropriate transfer func- 
tion (e.g. the radiation transfer function for the CMB 
temperature anisotropy), they can be used to constrain 
the value of the coupling parameters /nl and $nl- No 
significant detection of primordial non-Gaussianity has 
been reported from measurements of the three-point cor- 
relation function of the cosmic microwave background 
(CMB) anisotropies [HI, El El, El, El- The tightest 
limits are —4 < /nl < 80 at 95% confidence level 14 1. 
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If C(/nl) ~ O(snl) then the cubic correction should 
always be negligibly small compared to the quadratic 
one since curvature perturbations are typically 0(1O -5 ). 
However, this condition is not satisfied by some mul- 
tifield inflationary models such as the curvaton sce- 
nario, in which a large $nl and a small /nl can be si- 
multaneously produced. In this model, curvature per- 
turbations are generated by an additional scalar field, 
the curvaton, whose energy density is negligible dur- 
ing inflation [16|, E3, El EH- Non-Gaussianity is gener- 
ated by curvaton self-interactions which effectively con- 
tribute a non-quadratic term to the curvaton potential 
E2, [H [S3, UM- While the value and the sign of g NL 
depend upon the exact form of the self-interaction term 
(which can dominate the mass term if the curvaton mass 
is small enough and the curvaton vacuum expectation 
value during inflation is large enough [25j), it is generi- 
cally of magnitude |<?nl| ~ 10 4 — 10 5 for realistic models 
in which the ratio of the energy density of the curva- 
ton to the total energy density at time of decay is small. 
There are other realizations where one can have large $nl 
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and small /nl [2&I23- In ekpyrotic and cyclic models, 
/nl typically is of the order of a few tens while <7nl is 
of the order of a few thousand (28[. If /nl were small, 
then the imprint of non-Gaussianity would be detected 
only in four-point statistics such as the CMB trispectrum 
[29L f30l. r3lj] . Thus far, no observational limits have been 
set on <7nl by measuring the CMB trispectrum [§2, HH • 
Nevertheless, since the current bound |/nl| ^100 implies 
a relative contribution for the quadratic term of ~ 0.1 per 
cent, a third order coupling parameter |<7nl| ~ 10 6 should 
also be consistent with the data. 

Large-scale structures offer another venue to test for 
the presence of primordial non-Gaussianity. Deviation 
from Gaussianity can significantly affect the high mass 
end of the mass function [34], [351], the large-scale two- 
point correlation [3(| [13] , the bispectrum (38|, HiJ H3, E3 
of dark matter haloes hostin g th e observed galaxies as 
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42 , 43 1 and topological mea- 
451 ] . Recently, references 



well as void abundances 
sures of the cosmic web 

pel l47l |48| showed that the local quadratic coupling 
/nl^ 2 induces a scale-dependent bias Ab K (k, /nl) m the 
large-scale power spectrum of biased tracers, 



Ab R {kJ NL )=3f NL [b(M)-l]S c 



k 2 T(k)D(z) 



(2) 



where b(M) is the linear bias parameter, Hq is the Hub- 
ble parameter, T(k) is the matter transfer function nor- 
malised to unity as k — > 0, D(z) is the growth factor nor- 
malised to (l + z) -1 in the matter era and S c ~ 1.68 is the 
present-day (linear) critical density threshold. Reference 
[481 ] applied Eq. @ to constrain the value of /nl us- 
ing a compilation of large-scale structure data and found 
—29 < /nl < +69 at 95% confidence. These limits are 
comparable with those from the CMB, demonstrating the 
competitiveness of the method. Forthcoming all sky sur- 
veys should achieve constraints of the order of /nl ~ 1 
[ig . I49I IHol . [EH and should be sensitive to a possible scale- 
dependence of /nl [52]. On the numerical side however, 
while simulations of structure formation have confirmed 
the scaling Ab K (k, /nl) with k [H, [H, [13, HH , the exact 
amplitude of the non-Gaussian bias correction remains 
somewhat debatable. 

All numerical studies to date have only implemented 
the quadratic term /nl^ 2 - The purpose of this work is to 
quantify the impact of the cubic term ^nl^ 3 on the mass 
function and bias of dark matter haloes extracted from 
cosmological simulations and assess the ability of forth- 
coming measurements of the large-scale bias of galax- 
ies/quasars to constrain the size of a local cubic correc- 
tion. This paper is organized as follows. We begin with 
a brief description of the N-body simulations and illus- 
trate the extent to which the coupling (?nl0 3 affects the 
matter power spectrum and the halo mass function (Sec. 
IT!)) . We pursue with the non-Gaussian halo bias (Sec. 
IIIip . to which we derive analytically the scale-dependent 
and scale-independent contribution, Ab K and Abi, and 
demonstrates the large suppression of the simulated Ab K 
relative to theory. We then place limits on the coupling 



parameter <?nl and forecast constraints from future large- 
scale surveys and CMB experiments (Sec. IIVj) . We also 
show that our findings consistently apply to more gen- 
eral models with non-zero /nl and <?nl (Sec. [Vj) . We 
conclude with a discussion of the results in Sec. [VT] 



II. THE NON-GAUSSIAN SIMULATIONS 
A. Characteristics of the N-body runs 

We utilize a series of large N-body simulations of 
the ACDM cosmology seeded with Gaussian and non- 
Gaussian initial conditions. The (dimcnsionlcss) power 
spectrum of the Gaussian part <^>(x) of the Bardeen po- 
tential is the usual power-law A 2 ,{k) = k 3 P$(k) / (27r 2 ) = 
A < p(k/kQ) ns ^ 1 . The non-Gaussianity is of the "local" 
form $ = cp + .gNL</> 3 - We adopt the standard (CMB) 
convention in which ^(x) is primordial, and not extrap- 
olated to present epoch. It is important to note that the 
local transformation is performed before multiplication 
by the matter transfer function. T{k) is computed with 
CMBFAST [El] for the WMAP5 best-fitting parameters 
: h = 0.7, n m = 0.279, fl h = 0.0462, n s = 0.96 and 
a normalisation of the Gaussian curvature perturbations 
A, p = 7.96 x lO" 10 at the pivot point fc = 0.02MPC" 1 . 
This yields a density fluctuations amplitude a$ ~ 0.81 
when the initial conditions are Gaussian. Five sets of 
three 1024 3 simulations, each of which has <7nt, = 0, ±10 6 , 
were run with the N-body code GADGET2 [53] . We used 
the same Gaussian random seed field cj> in each set of 
runs so as to minimise the sampling variance. We also 
explored scenarios with non-zero /nl and c/nl and ran 2 
realisations for each of the non-Gaussian models charac- 
terized by (/nlj <?nl) = (±100,-3xl0 5 ). In all cases, the 
box size is 1600 /i _1 Mpc with a force resolution of 0.04 
times the mean interparticle distance. The particle mass 
of these simulations thus is 3.0 x 10 11 M©//i, enough to 
resolve haloes down to 10 13 M Q //i. 

In the curvaton scenario, generic polynomial interac- 
tion terms of the form Xm^{u /m (7 ) n (where A is some 
coupling strength and m a is the curvaton mass) to the 
quadratic potential of the curvaton field a yield |<?nl| ^ 1 
even when the non-linearity parameter /nl is very small 
[HI [12. One typically finds |g^ L | - O(10 4 ) - O(10 5 ) 
when /nl varies in the range —100 < /nl < 100. For 
practical reasons however, the values of <;nl adopted in 
our simulations are about an order of magnitude larger 
so as to produce an effect strong enough to be unam- 
biguously measured despite the small simulated volume. 
Furthermore, we have also considered positive and neg- 
ative values of <?nl so as to assess the sensitivity of the 
non-Gaussian bias to the sign of the coupling parameter. 
The simulations with (/nl, Snl) = (—100, —3 x 10 5 ) may 
be seen as a particular realisation of the curvaton model 
in which the coupling constant A is positive, and the non- 
quadratic term is very steep (n ~ 5 — 10) but contributes 
little to the total curvaton potential. 
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B. Properties of the initial density field 

In order to ensure that the initial conditions are suc- 
cessfully generated, we measure at the redshift of our 
initial conditions, z = 99, the (normalized) skewness 
S 3 (R,z) = (S 3 R z )/a 4 and kurtosis S 4 (R,z) = ((S 4 R z )/ - 
3a 4 ) /a 6 of the density field 5r. z smoothed with a (spheri- 
cally symmetric) window function of characteristic radius 
R. We adopt a tophat filter throughout this paper. Note 
also that a(R, z) is the variance of smoothed density fluc- 
tuations at redshift z. 

In the weakly nonlinear regime, the skewness and kur- 
tosis of the density field may be written as the sum of 
a part generated by gravitational clustering and a part 
induced by primordial non-Gaussianity. For Zel'dovich 
initial conditions [58j and Q m (z) « 1, the contribu- 
tion generated by gravitational instabilities reads as 

M M M M 



f el (i?, z) = 4 - (n cff + 3) (3) 

, n , 272 50 , oN 7 . „, 2 

4 (R,z) = — — (n e ff + 3) + - (n e ff + 3) , 



S 



5 Zcl 




R / h-'Mpc 



where n c g(R) is the effective spectral index at the 
smoothing scale R, 



n cS (R) 



d\n<j 2 (R, z) 
dlni? 



(4) 



Note that the initial skewness and kurtosis given by the 
Zel'dovich approximation differs from the exact values 
predicted by perturbation theory, to which they asymp- 
tote in the limit D(z) — *■ oo 63). In addition, the cubic 



coupling <?nl0 3 induces a nonzero kurtosis Sf ri (R,z) 



FIG. 1: Skewness and kurtosis of the initial (z — 99) density 
field as a function of smoothing radius. While the top panel 
shows the sum of the contributions arising from the Zel'dovich 
dynamics and from primordial non-Gaussianity, S3 = S3 01 + 
S^ T \ the bottom panel only shows the absolute value of the 
primordial kurtosis, |S4 ri |. Symbols represent the numerical 
results averaged over the realisations. They have been slightly 
shifted horizontally for clarity. Error bars show the scatter 
among the realisations for the models with (Jnl = 10 6 . Solid 
lines indicate the theoretical expectations. 



9NhS^\R,z) at leading order which can be computed absolute difference |5f ri | between the kurtosis in the non- 



analytically from the relation 



a^Pi^z) 




where 



a R (k,z) 



3r^ m i?o 



•kn-i|,z) , (5) 



D(z)k 2 T(k)W R {k) (6) 



is evaluated at redshift z and W R (k) is the Fourier trans- 
form of the tophat function. When /nl = 0, primordial 
skewness is not generated at the first order and, there- 
fore, may be neglected. 

Fig. [T] displays the initial skewness (top panel) and kur- 
tosis (bottom panel) obtained upon distributing the dark 
matter particles onto a regular 512 3 mesh (i.e. of cell size 
w 3 /i _1 Mpc) using the cloud- in-cell (CIC) interpolation 
scheme. Symbols represent the numerical results aver- 
aged over the realisations. Because Sf cl (which is the sole 
contribution to the kurtosis in the Gaussian case) varies 
considerably between the realisations, we only show the 



Gaussian (^nl = ±10 6 ) and the Gaussian (#nl = 0) 
runs. As expected, the skewness is very similar among 
the Gaussian and non-Gaussian simulations. While S3 
in the simulations agrees well with the skewness induced 
by the Zel'dovich dynamics Eq. §5§ (solid curve), it is 
gradually suppressed as the filtering radius R approaches 
the cell size, presumably because of the finite resolution 
which smoothes also the fluctuations. Note that, in the 
simulations with nonzero <7nl and /nl, there is a large 
primordial skewness which is of magnitude | Sf rl | ~ a few 
on scale R < 50 /i~ 1 Mpc, consistent with theory. As can 
also be seen, the absolute value of the primordial kurtosis 
increases sharply with R, in very good agreement with 
the theoretical prediction. 



C. The matter power spectrum 

Non-Gaussian corrections to the primordial curvature 
perturbation can renormalise the input power spectrum 
of fluctuations used to seed the simulations. Since 
our simulations implement its unrenormalised version 
A?(fc) = ^(fc/fco)™^ 1 , it is desirable to ascertain the 
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effect of the local coupling term on the simulated density 
power spectrum before discussing the halo mass function 
and bias. For /nl models with |/nl| ^ 100, renormal- 
isation effects are unlikely to be noticeable due to the 
limited dynamical range of current cosmological simula- 
tions [64]. As we will see shortly however, they can be 
significant in simulations of (/nl models with similar level 
of non-Gaussianity. 

The cubic order term ^nl^ 3 renormalises the ampli- 
tude A$ of the power spectrum of initial curvature per- 
turbations to — > + 6(7nl(0 2 ), where 



d 3 k 



(7) 



For scale invariant initial conditions, ((f> 2 ) has a logarith- 
mic divergence at large and small scales (see [64| for a 
more detailed discussion of perturbative corrections in 
non-Gaussian cosmologies). In practice, a low- and high- 
k cutoff are naturally provided by the finite box size and 
the resolution of the simulations. Therefore, the effec- 
tive amplitude of density fluctuations in non-Gaussian 
simulations with cubic coupling is erg + Sas with 



5a & = 3g N L(^ 2 ) 



(8) 
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k ■ 



1— n. 
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l-n s 
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Recall that ko = 0.02 /iMpc -1 is our choice of normalisa- 
tion point, and fc m ; n and fc max are the integration limits 
set by the fundamental mode and the Nyquist frequency 
of the periodic cubical box over which the initial condi- 
tions are generated. Equivalently, 



<Ws = 3<7nl ( J [l 



N n s 



A r! 



1 



(9) 



where N = 1024 is the number of mesh points along one 
dimension. This result becomes 5a% = 3(?nl ^0-(N)A ( p in 
the scale-invariant limit n s — > 1. For the cosmological 
setup considered here, the absolute deviation is 



<5cr 8 « 0.015 



(9mA 
V 10 6 / 



(10) 



This correction is fairly large for the values of <7nl 
adopted here and, therefore, must be taken into account 
in the comparison between the theory and the simula- 
tions. As we will see below, this is especially important 
when studying the high mass tail of the halo mass func- 
tion which is exponentially sensitive to the amplitude of 
density fluctuations. 

The cubic coupling term c/nl^ 3 can also induce a 
scale-dependent correction to the matter power spec- 
trum which can be quantified by the fractional change 

/?m(fc,3NL) = Pmm(k, #Nl) / -P m m(fc, #NL = 0) - 1. In 

Fig. [H symbols show the result of measuring /3 m (fc, (/nl) 
from the snapshots at z = and 2 after correction of 
the normalisation shift \28<j%/<j%\ — 0.037. There is some 



evidence for a scale-dependent correction at wavenumber 
k > 0.1 /iMpc -1 but the resulting deviation is broadly 
consistent with zero. We will thus neglect /3 m (fc,(/NL) 
henceforth. 



D. The halo multiplicity function 

Haloes were identified using the MPI parallelised ver- 
sion of the AHF halo finder [65| which is based on the 
spherical overdensity (SO) finder developed by [fjfj]. The 
main reason for using a SO finder is that it is more 
closely connected to the predictions of the spherical col- 
lapse model, on which most of the analytic formulae pre- 
sented in this paper are based. Namely, the virial mass 
M of a halo is defined by the radius at which the inner 
overdensity exceeds A v ; r (z) times the background den- 
sity p(z) [67], [Hj]. The value of the overdensity thresh- 
old A v ; r (z) is obtained from the collapse of a spherical 
tophat perturbation and has a dependence on redshift 
through the matter density £l m (z) [iil, [7(|. We discard 
poorly resolved haloes and only study those containing 
at least 34 particles or, equivalently, with a mass larger 
than M = 10 13 M Q /h. 

A naly tic arguments based on the Press-Schechter the- 
ory [Tl|, l72J predict that the halo mass function n(M, z) 
is entirely specified by the distribution v]{y) of first- 
crossings, or multiplicity function 



vf(v) = M 



2 n(M,z) d\nM 



(11) 



The peak height v(M,z) — S c (z)/a(M) is the typical 
amplitude of fluctuations that produce haloes of mass 
M by redshift z. Here and henceforth, <j(M) denotes 
the variance of the density field 5m smoothed on mass 
scale M oc R 3 and linearly extrapolated to present epoch, 
whereas 5 c (z) « 1.68D(0)/ D(z) is the critical linear over- 
density for (spherical) collapse at redshift z. 

Despite the lack of a compelling theoretical descrip- 
tion of the multiplicity function for Gaussian initial con- 
ditions, the fractional deviation from Gaussianity can be 
modelled accurately using the Press-Schechter formalism. 
In this approach, the halo mass function n(M, z) is re- 
lated to the probability P{> 5 C , M) that a region of mass 
M exceeds the critical density for collapse 5 c (z) through 
the relation n(M,z) = -2 (p/M) dP/dM. The non- 
Gaussian fractional correction to the multiplicity func- 
tion then is R(v,g NL ) = f{v,gm)/f{y,0) = (dP/dM)(> 
5 c ,M,g NL )/(dP/dM)(> 6 C ,M,0). The level excursion 
probability P(> 5 C ,M, <?nl) can be computed once the 
probability distribution function (PDF) of the smoothed 
density field 5m, P{5m), is known. Here, we will consider 
the simple extensions proposed by [73[ and [jH , in which 
P{5m) is generically expressed as the inverse transform 
of a cumulant generating function. Both extensions have 
been shown to give reasonable agreement with numerical 
simulations of non-Gaussian cosmologies [53l l55l. |75| . 
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FIG. 2: Non-Gaussian fractional correction /3 m (&, <7nl) = 
Pmm(k, g^h)/ Pmm(k,0) — 1 to the matter power spectrum 
after subtracting a scale-independent normalisation shift 



"} induced by the cubic coupling (/nl</ 



In [73| , the saddle-point technique is applied directly to 
P(Sm)- The resulting Edgeworth expansion is then used 
to obtain P(> <5 C , M, <?nl)- For /nl non-Gaussianity, ref- 
erence 53] found that the resulting non-Gaussian mass 
function agrees well with the simulations. For <7nl non- 
Gaussianity, neglecting cumulants other than the kurto- 
sis S^ n (M) (Hereafter, we drop the superscript for con- 
ciseness) and truncating the series expansion at Si(M), 
the non-Gaussian fractional correction reads 



«{l + i^S 4 (^-to 2 -3) (12) 
= jl + ^CT 2 S 4 (^ 4 -6^ 2 + 3) (13) 



{v 2 - 3) \ exp [u 2 5a 8 ] 



1 d{a 2 S 4 ) 
4! d\nv 



after integration over regions above 8 c (z). Note that we 
have omitted writing the redshift dependence explicitly. 
Strictly speaking however, R(v, (?nl) depends distinctly 
upon the variables M (or v) and z due to the presence of 
a 2 S4(M). Our notation is motivated by the fact that the 
measured non-Gaussian correction, as plotted in Figf?) 
appears to depend mostly on the peak height. The ex- 
ponential factor in the right-hand side is the correction 
induced by the renormalisation of the amplitude of lin- 
ear density fluctuations, Eq. (Ti"Tj|) . For consistency, we 
have also used the Press-Schechter multiplicity function 
to derive this last term although a Sheth-Tormen mass 
function [76[ may be more appropriate. 

In [74j], it is the level excursion probability P(> 5 C , M) 
that is calculated within the saddle-point approximation. 
Including only a cubic coupling 5nl0 3 and truncating the 



resulting expression at the kurtosis, we find 



P(> (5c,M,<?nl) 
1 a 



(14) 



V27T S, 

x exp 



1 + 3#nlo- 



2/J.2\ 



2(7 2 



I-65 



NL \ 



^ 4 ^2 r2 

12 

' 12 c 



at first order in g^L ■ Note that we have already included 
the renormalisation of the fluctuation amplitude. For 
rare events, it«1 and the first parenthesis in the right- 
hand side can be neglected. To ensure that the resulting 
mass function is properly normalised, we follow [77| and 
use 



!/*/(!/*) = M 



2 n(M, z, <?nl) din M 



din v+ 



(15) 



for the non-Gaussian mass function, where = <5*/er, 
S+ = S c y/l — 2<5ct8 — 5*4(52/12 1 and / is the same mul- 
tiplicity function as in the Gaussian case. Taking the 
derivative of the level excursion probability then gives 



(dP/dM){> <5 c ,M,g NL ) 
(dP/dM)(> 5 C ,M,0) 



exp 



4!^ + V S<78 
5* 1 S 3 r dS 4 

JL _| S. i_ 

S c 4! (5* diner 



(16) 



The fractional change in the multiplicity function even- 
tually reads as 



-Rmvj^Snl) « exp 

xfl 



4! 



a S4 + v Sag 



(17) 



-a 2 S A 



v 2 d(a 2 S 4 ) 
4! d\ni> 



after expanding 5* at the first order and ignoring the shift 
in the normalisation amplitude, i.e. 6+ 1 w 1 =p S*4(5 2 /24. 
In the limit <7 2 5 4 <C 1 and v ^S> 1, the two theoretical ex- 
pectations reduce to 1 + ^ 4 cr 2 5'4/24. However, they differ 
in the coefficient of the v 2 a 2 Si term, which is —1/4 and 
— 1/8 for the LV and MVJ formula, respectively. Thcrc- 



1 In local /nl models, this formula only involves the skewness 
53. As it is incorrectly quoted in some of the literature on non- 
Gaussian halo mass functions, let us write down its explicit ex- 
pression: 

52 dS 3 i 



i?MVj(^, /nl) = exp 



S 3 8 3 c 



6cr 2 



6(5* din a 



or, in terms of the peak height v = S c /cr 
,3 



-Rmvj(^> /nl) ~ exp 



-0-53 



1 



V v d(aS 3 ) 

— (703 

3 6 d lnf 



after expanding 5* = ^/l — S 3 8 c /3 at the first order. 
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FIG. 3: Variance a (dotted), skewness aS 3 (dashed) and 
kurtosis a 2 S^ (solid) of the smoothed linear density field 
8m as a function of mass scale M. 



fore, we shall also consider the approximation 



R(v,gNh) = cxp 



— a 2 Si + v 2 &Of, 



(18) 



1 - -^<T 2 S 4 - 



v 2 d{a 2 Si) 
4! d\nv 



which is designed to match better the Edgeworth expan- 
sion of [7|| when the peak height is v ~ 1. 

Calculating the fractional change in the mass function 
requires knowledge of the kurtosis Sa(M) = 
of the smoothed density field 5m , which we compute an- 
alytically using the general formula (valid for n > 3) 
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2 3 4 5 6 7 

FIG. 4: Top panel : Fractional correction to the Gaussian mul- 
tiplicity function of dark matter haloes as a function of the 
peak height v(M,z) for a coupling parameter £nl = ±10 6 . 
The dotted, dashed and solid curves show the theoretical pre- 
dictions eqs. (I13p . (|17[1 and (|18[1 at z — 0, respectively. Error 
bars denote Poisson errors. For illustration, M — 10 15 Mq/A 
corresponds to v = 3.2, 5.2, 7.7 at redshift z = 0, 1 and 
2, respectively. Similarly, M = 10 14 M Q /ft and 10 13 M Q /ft 
correspond to v — 1.9, 3, 4.5 and 1.2, 1.9, 2.9 respectively. 
Bottom panel : a comparison with equation (|18p evaluated at 
z = and 2. 



a 2 - 2 s£\M) = n! (n/^«*r(**)^(*«)l 

x a M (|ki + ---+k n _i |) , (19) 

where a.M{k) = &n(k, z = 0) in what follows. Over the 
mass range probed by our simulations, 10 13 < M < 
5 x 10 15 M Q /h, the normalised kurtosis cr 2 ^ 1 ^/) is 
a monotonic decreasing function of M that varies in the 
narrow range ~ 4— 6x 10 -7 for the top-hat filter assumed 
here (see Fig. [3]). Note also that the <j 2 S± term dominates 
the total contribution to the non-Gaussian correction eqs. 
(fl~3|) . p7| and (JTHJ) when the peak height is v > 2 (One 
finds |d(o- 2 ,5*4)/dln^| < O.Ict 2 ^). 

The fractional correction is plotted in Fig. 2] for the 
haloes extracted from the simulations at redshift z = 0.3, 
0.5, 1 and 2. In the top panel, the data are compared 
to the theoretical predictions eqs. (fTB"]) . (fTT]) and (fl8|) 
evaluated at z = 0. As we can see, the level of non- 



Gaussianity in the halo multiplicity function is consis- 
tent with the theory. Our approximation (|18p performs 
equally well regardless of the sign of #nl ■ It agrees better 
with the measurements than the formulae of [74| which 
significantly overestimates the data for <?nl = 10 6 , and 
than that of [73| which is not always positive definite for 
<7nl = — 10 6 . The bottom panel shows that the discrep- 
ancy somewhat worsens at higher redshift, especially in 
the case <7nl = 10 6 . However, it is possible the agree- 
ment may be improved by adding higher order powers of 
a Si and higher order cumulants. 



To conclude this section, one should keep in mind that 
all these extensions are based on Press-Schechter and, 
therefore, provide a bad fit to the Gaussian mass function 
of haloes. In this respect, excursion set approaches may 
be more promising since they seem to reproduce both the 
Gaussian halo counts and the dependence on /nl [zl, Isl • 
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III. THE NON-GAUSSIAN BIAS SHIFT 



Theoretical considerations 



In order to calculate the scale-dependent bias correc- 
tion induced by the <7nl coupling term to the correlation 
of haloes of mass M collapsing at redshift z, we follow 
p7l ] and consider the two-point correlation £hh(r) of re- 
gions of the smoothed density field 8m above a threshold 
S c (z) = i>{z)a. The two-point correlation function of this 
level excursion set, which was originally derived by [37] |. 
can be expressed in the high threshold approximation as 



00 n— 1 



^-^ j\(n — j)\ 

n=2 j=l J y J 1 
(n)( Xi,--- ,Xl, X 2 ,--- ,X 2 

j times (n — j) times 



(20) 



where r = X! — x 2 . For the non-Gaussian model consid- 
ered here, the leading-order correction induced by non- 
zero three-point and four-point correlations of the density 
field reads 



A6 a h = -^ (3) (xi,x 1 ,x 2 ) + ^ i 

+ i^ (4) ( Xl ' Xl ' X2 ' X2 ) 



t:£ (4) (xi,Xl,Xl,X 2 ) 



(21) 



The non-Gaussian correction APhh to the power spec- 
trum of biased tracers is obtained by Fourier transform- 
ing this expression. 

In the case /nl = and c?nl 7^ 0, only the four- 
point functions contribute at first order. It should also 
be noted that, at linear order, £( 2 )(xi,x 2 ) amounts to 
a renormalisation of the linear bias and, therefore, does 
not contribute to the scale-dependent correction. Details 
of the calculation can be found in Appendix lAl In short, 
the non-Gaussian correction APhh in the limit of long 
wavelength k <C 1 is given by the Fourier transform of 
i^ 4 £ (4) (xi , xi, x 2 , x 2 )/3cr 4 , 

AP hh (fc) = 5N iV S^\M)a M {k)P^k) (22) 
= g Nh b 2 (z)6 2 c (z) 5< 1) (M)a M (fc)P (fc) , 



where we have used b^(z) — v 2 /S c (z) as is appropriate 
for high density peaks. The smoothing window that ap- 
pears in «A/(fc) effectively makes little difference because 
we are considering the limit where fc _1 is much larger 
than the smoothing radius, so we will omit it in the fol- 
lowing. For small non-Gaussianity, we can also write 
AP hh ~ 2b L Ab K P s {k) where P s (k) = a^P^fc) is the 
power spectrum of the smoothed density field. The scale- 
dependent bias correction Ab K (k, <?nl) can eventually be 



z=0.51 
b(M) = 2.52 




0.005 0.01 



0.05 0.1 



k / hMpc" 



FIG. 5: Top panel : Non-Gaussian bias correction com- 
puted from the halo-matter power spectrum of haloes of 
mass M > 2 x 10 13 Mq//i extracted from the snapshot 
at z — 0.5 (filled symbols). The solid curve represents 
Pmh(k, gNh)/ Pmh(k, 0) — 1 with a non-Gaussian bias shift 
Ab(k, <?nl) given by Eq. ([26}. The dashed, dotted and dotted- 
dashed curves show three separate contributions that arise at 
first order in <7nl- Bottom panel : A&(&,<?nl) is replaced by 
the theoretical model Eq. ((27}. The shaded region indicates 
the data points used to fit the parameters e K and e\. Error 
bars indicate the scatter among 5 realisations. 



recast into the form 



Ab K (k,g NL ) = -g NL b L (z)5 2 (z) sP(M)a^(k) (23) 



- A 9^5 c {z) SP(M) Ab K (k, /nl - 1) , 



where Ab K (k, /nl) is the scale-independent bias induced 
by the quadratic coupling /nl^ 2 , Eq. (J2J). Wc have 
also assumed the Eulerian bias prescription b(M) = 
l + b L {M). 

The change in the mean number density of haloes also 
creates a scale-independent shift which we denote by 
A&i(<7nl)- As shown in [{53| for /nl models, the inclusion 
of this correction noticeably improves the agreement with 
the simulations at wavenumber k < 0.1 /iMpc -1 . Using a 
peak-background split and considering the limit of small 



8 



non-Gaussianity, this contribution reads 
1 d 

A&i(£Tnl) = ^-ln[i?(^,ffNL)] 



4! v ' d\sxv 



3i/) ct 2 ^ + 2^cr 8 (24) 

d(<T 2 S 4 ) V d 2 (<T 2 S A )~ 



4! din; 



after truncating A&i at first order in <7nl This approxima- 
tion should perform reasonably well for moderate values 
of the peak height, v < 4, for which the fractional change 
in the mass function, equation (|18p . matches well the nu- 
merical data. It is worth noticing that A&i(<?nl) has a 
sign opposite to that of <?nl because the bias decreases 
when the mass function goes up. A&i(<7nl) also includes 
a correction induced by the renormalisation of erg. In 
practice, we estimate A&i(<7nl) for a given halo sample 
by evaluating a 2 S^ and v at the scale corresponding to 
the average halo mass M of the sample, as it is unclear to 
which extent R(v, <?nl) agrees with the data in the limit 
M > M*. 



<3 



0.5 - 



0.5 



- 



0.5 

-1 



a ' lO'^M^SxlO^h'-'Mg,' gl!=-10 6 
* M>2xl0 13 h-'M S , g NL =-10 6 



'I I I" 



i i i 'i 



I I I I I I I I I III II II III llll I 



M>2xl0 13 h-'M G , g NL =10 6 
10 13 <M<2xl0 13 h- 1 M o , g NL =10 6 



il I In 



3 

b(M) 



6 7 



B. Comparison with the simulations 

To assess the effect of primordial non-Gaussianity on 
the halo bias, we will consider the ratios 



-Pmh(fc,ffNL) 

P m h(M) 
-Phh(fc,ffNL) 

P hh (fc,0) 



1 



Afr(fc,ff NL ) 

b(M) 



.<5cr 8 



CT8 



(25) 



Ab(fc,g N L) 
b(M) 



where A6(fc,<?NL) is generally the sum of a scale- 
dependent and a scale-independent term. One should 
bear in mind that the scale-independent shift 2Sas/ag 
arises from the matter power spectrum and, therefore, is 
distinct from the term —Ivba^ja appearing in A6i. Fol- 
lowing [55j | , we shall also quantify the departure from the 
theoretical scaling as a function of wavemode amplitude 
with the ratio A6 S /A6', where A6 S is the non-Gaussian 
bias correction measured from the simulation and A6* is 
Eq. 

We interpolate the dark matter particles and halo cen- 
tres onto a regular cubical mesh. The resulting dark 
matter and halo fluctuation fields are then Fourier trans- 
formed to yield the matter-matter, halo-matter and halo- 
halo power spectra P mm (k), -fmh(fc) and Phh(fc), respec- 
tively. These power spectra are measured for a range of 
halo masses and redshifts, covering the relevant range of 
statistical properties corresponding to the available data 
sets of galaxy or quasar populations with different lumi- 
nosities and bias. Note that these quantities are com- 
puted on a 512 3 grid, whose Nyquist wavenumber is suf- 
ficiently large (rj 1 /iMpc -1 ) to allow for an accurate 
measurement of the power in wavemodes of amplitude 
k < 0.1 /iMpc -1 . The halo power spectrum is corrected 



FIG. 6: Best-fitting e K and Abi + e\ as a function of halo bias 
and gfNL = ±10 6 , for two different mass cuts as indicated in 
the Figure. In the top panel, the dotted curve is our best fit 
to e K (6, <?nl). In the bottom panel, the solid and dashed lines 
show the scale-independent shift Abi predicted by a peak- 
background split, Eq. (|24|l . for the low and high mass samples, 
respectively. 



for the shot-noise due to the discrete nature of dark mat- 
ter haloes, which we assume to be the standard Pois- 
son term 1/nh- This discreteness correction is negligi- 
ble for P mm (k) due to the large number of dark matter 
particles. Yet another important quantity is the linear 
halo bias b(M) which must be measured accurately from 
the Gaussian simulations as it controls the magnitude of 
the scale-dependent shift. Here, we shall use the ratio 
Pmh(k) / P m m(k) as a proxy for the halo bias since it is 
less sensitive to shot-noise. 



1. An effective non-Gaussian bias correction 

Summarizing the analytical considerations of Sec. 
IIII A[ non-Gaussianity of the <7nl type add a correction 
A6(/c,(7nl) to the bias b(k) of dark matter haloes which 
is at leading order 



Ab{k,g NL ) = A6 K (fc,g N L) + A6i(g N L) 



(26) 



We found, however, that this theoretical expectation 
significantly overestimates the magnitude of the non- 
Gaussian bias shift measured from the simulations. 
This is exemplified in the top panel of Fig. O where 
Pmh(k,gNh)/Pmh(k,0) — 1 is plotted for haloes of mass 
M > 2 x 10 13 Mo/h identified at z = 0.5. Clearly, the 



9 




FIG. 7: Non-Gaussian bias correction measured in the simulation outputs at redshift < z < 2 for haloes of mass M > 
2 x 10 13 MQ/h. In each panel, the upper plot shows the ratio Phh(k, gwh) / Phh(k, 0) — 1 (dashed curves, empty symbols) 
and P m h(fc, <7NL)/Pmh(fc, 0) — 1 (solid curves, filled symbols). The error bars represent the scatter among 5 realisations. The 
respective output redshift and linear halo bias are also quoted. The bottom of each panel displays the departure from the 
theoretical prediction, A& s /A&*. The shaded area indicates the domain where the deviation is less than 20 per cent. The 
parameters e K and ei are fitted individually to each sample. For illustration, e K takes the best-fit values 0.06 and 0.60 for the 
lowest and highest biased samples, respectively. 
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FIG. 8: Same as Fig. but for haloes in the mass range 10 13 < M < 2 x 10 13 M Q /h. 



predicted scale-dependent correction Ab K is much steeper 
than measured from the halo samples. In order to im- 
prove the agreement with the numerical data, we modify 
the above relation as follows : 

Ab(k, <?nl) = e K Ab K {k, g Nh ) + [A6i( 5nL ) + d] , (27) 

and treat e K and ei as free parameters that we fit to our 
measurements of the cross-power spectrum (weighted by 
the scatter among 5 realisations) in the range 0.005 < 
k < 0.03 /iMpc -1 where the scale-dependent effect is 
largest. The bottom panel of Fig. [5] shows the result- 
ing best-fit contributions e K Ab K and Ab\ + c\ for the 
halo sample mentioned above. As seen in Fig. e K and 
ei appear to depend mainly upon the linear halo bias 
b(M) and the coupling parameter <7nl, although depen- 
dencies on redshift or other halo observables are not ex- 
cluded (The data is too noisy for a reliable estimate of 
these). The most striking feature of Fig. [5] is the func- 
tional dependence of e K on b(M) and <7nl- Firstly, e K 
is a monotonically increasing function of the bias and 
never reaches unity, even for the most biased samples 
for which the high peak approximation should be valid. 



Secondly, e K is noticeably larger for <?nl = — 10 6 , suggest- 
ing thereby that second (and higher) order contributions 
to the scale-dependent bias may be important. Further- 
more, for b < 1.5 where the high peak approximation 
breaks down, there is some evidence that the effect re- 
verses sign. The bottom panel of Fig. [6] shows that the 
scale-independent correction has sign opposite to that of 
<?nl ; m agreement with theoretical expectations from the 
peak-background split (see Sec lIII A| ). However, whereas 
for b < 3 the magnitude of the correction is comparable 
to that predicted by Eq. (HM)) , it is considerably larger 
for b > 3, reaching up to 5-10 per cent of the linear halo 
bias. 

Assuming e K is a function of b(M) and <?nl only and 
asymptotes to a constant in the highly biased limit, we 
find that the following parametrised form 

e K {b, 5nl) = ci - c 2 ,9nl - — ^—rr (28) 

1 + CXp(C40) 

captures reasonably well the increase of e K with halo bias 
for 1.5 < b(M) < 7. The best-fit values of the parameters 
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are 

Cl = 0.59 ±0.03, c 2 = (6.0 ±0.9) x 10~ 8 (29) 
c 3 = 2.1±0.5, c 4 = 0.88 ±0.13 . 

We do not provide a fitting formula for e\ (or A&i + £i) 
since it is not directly measurable in real data. 

Since both the kurtosis of the initial density field and 
the mass function of dark matter haloes are in fairly good 
agreement with theoretical expectations, the discrepancy 
in the scale-dependent bias A& K (fc,<?NL) indicates that it 
is the high peak approximation considered here which is 
inaccurate, even in the limit b(M) ^> 1 where it is sup- 
posed to work best. Our perturbative approach may be 
one of the reasons for this mismatch. Namely, we have 
derived the scale-dependent bias A& K (fc,gNL) under the 
assumption that the non-Gaussian correction to the halo 
power spectrum is small whereas, for most of the bias 
range covered by our halo catalogues, the effect is in fact 
already very large at k < 0.01 /iMpc^ 1 . Note, however, 
that for the quadratic coupling /nl0 2 , this perturbative 
treatment predicts an effect of the right magnitude 
Explaining these findings clearly requires a better the- 
oretical understanding, which we leave for a future in- 
vestigation. Notwithstanding this, the phenomenological 
prescription (|2T[) with parameters e K and t\ fitted to the 
data provides, as we will see below, a good description 
of the large-scale halo power spectrum in simulations of 
5nl models. 



2. Non-Gaussian bias from auto- and cross-power spectra 

We have measured auto- and cross-power spectra for 
a range of halo masses and redshifts spanning the range 
< z < 2. The ratios defined in Eq. (|2"5)) are shown 
in Figs [7] and H as a function of wavenumber for the 
mass threshold M > 2 x 10 13 M Q /h and the mass bin 
10 13 < M < 2 x 10 13 M©//i, respectively. The fractional 
deviation A& s /A&* is also shown at the bottom of each 
panel. The shaded region indicates a departure less than 
20 per cent. Error bars denote the scatter around the 
mean and, therefore, may underestimate the true errors 
as they are computed from a small number of realisations. 
Note that, in order to reduce the impact of sampling vari- 
ance, we first compute the ratios P m h(fc, <7nl)/ Pmh{k,0) 
and Phh{k, gNh)/Phh(k, 0) for each realisation before cal- 
culating the average. 

As we can see, once e K and e\ are fitted to the ratio of 
cross-power spectra, the theoretical prediction Eq. (j2"?| 
provides a reasonable description of the non-Gaussian 
bias in the halo power spectrum Phh, indicating that 
non-Gaussianity does not generate much stochasticity 
and the predicted scaling Ab K (k,g^) oc k~ 2 T{k)~ 1 ap- 
plies equally well for the auto- and cross-power spectrum. 
This was also found to be true in / NL models [53]. The 
inclusion of a scale-independent correction iSb\ + ei sig- 
nificantly improve the agreement at k < 0.03 /iMpc -1 . 



For the highly biased samples b > 4 however, this cor- 
rection is so large that the non-Gaussian bias shift re- 
verses sign at wavenumber k > 0.05 hMpc^ 1 . Such 
an effect is not seen in the simulations, but we expect 
large deviations from the relation (|27[) in that range 
of wavenumber, where second- and higher-order correc- 
tions induced by the cubic coupling .9nl</> 3 together with 
the nonlinear bias created by the gravitational evolu- 
tion of matter density fluctuations may become impor- 
tant. Even though the data is noisier due to the low 
number density of haloes, it is worth noting that, for 
the highly biased samples at k < 0.01 /iMpc -1 , the 
cross-power spectrum F m h(fc,.9NL = — 10 6 ) goes neg- 
ative while -Phh(fc,<?NL = — 10 6 ) remains positive and 
increases sharply, in agreement with the analytic pre- 
diction. Still, there is some evidence that the ratio 
Phh(k,gNL = — 10 6 )/Phh(fc,0) — 1 saturates at a value no- 
ticeably larger than -1 before the sharp upturn, whereas 
our model predicts 28(7%/ cr% — 1 « —0.96 at the minimum. 

Fig. [5] further explore the effect in the low mass sam- 
ples, for which the z = haloes with b(M) » 1.15 consti- 
tute an almost unbiased sample of the density field. In 
this case, the sign of the scale-dependent contribution is 
reversed, namely, the large-scale halo power spectrum in 
simulations of c/nl = — 10 6 is enhanced relative to that 
of the Gaussian ones. This is in rough agreement with 
the theory, which predicts a similar effect for b(M) < 1. 
Again, haloes with a similar bias also have a compara- 
ble scale-dependent and scale-independent bias regard- 
less of mass or redshift. Finally, note that the sample at 
z = 1.39 shown in Fig.[5]corresponds closely to the quasar 
sample used by [48|], for which z = 1.8 and b = 2.7. 

IV. CONSTRAINTS ON THE COUPLING 
PARAMETER #nl 

A. Constraints on <?nl from current large-scale 
structure data 

Reference 48] took advantage of the scale-dependence 
of the bias to constrain /nl from a sample of highly bi- 
ased luminous red galaxies (LRGs) and quasars (QSOs). 
It is straightforward to translate their 2-a limit —29 < 
/nl < +69 into a constraint on <?nl since the non- 
Gaussian scale-dependent bias A6 K (fc, (?nl) has the same 
functional form as Ab K (k, /nl)- 

Constraints will arise mostly from the QSO sample at 
median redshift z = 1.8, which covers a large comoving 
volume and is highly biased, b = 2.7. In light of our 
results (see Fig. EJ), we expect the parameter e K (b, <7nl) 
to vary with <?nl- However, in order to simplify the 
analysis, we will assume that, at fixed b, e K (b, 3nl) is 
given by the mean of e K (b, #nl = ±10 6 ). For a sam- 
ple with bias b ~ 2.7, this implies e K ~ 0.4. Further- 
more, assuming the typical mass of QSO-hosting haloes 
is ~ 10 13 Mq/A yields s£\m) ~ 2.3 x lO" 4 . Hence, 
the multiplicative factor (1/4) (5 c (z)e K S , 3 1 \M) is approx- 
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imately ~ 8.4 x 10 5 . Our limits on c/nl thus are 

- 3.5 x 10 5 < 5 NL < +8.2 x 10 5 (30) 

at 95% confidence level. The scale-independent correc- 
tion A&i + ei is not directly measured as it adds up to 
the bias b which is fitted to the data. For the limits ob- 
tained here, \Ab\ + ei| should be much smaller than b and 
can thus be ignored. Note also that, whereas the non- 
Gaussian bias scales as D{z)~ 1 in / NL models, we have 
Ab(k, <?nl) fx D(z)~ 2 for <7nl non-Gaussianity, so one can 
achieve relatively larger gains from measurements of high 
redshift tracers. In fact, the extent to which one can im- 
prove the observational bounds will strongly depend on 
our ability to minimize the impact of sampling variance 
caused by the random nature of the wavemodes, and the 
shot-noise caused by the discrete nature of the tracers. 
By comparing differently biased tracers of the same sur- 
veyed volume [5l| and suitably weigh ting galaxies (e.g. 
by the mass of their host halo) [80, l81j |. it should be 
possible to circumvent these problems and considerably 
improve the detection level. 



B. Predictions for future LSS surveys 

References [5(1 HH, H3, [HJ applied the Fisher matrix 
formalism to forecast constraints on /nl from forthcom- 
ing galaxy redshift surveys. Here, we will simply try to 
estimate the detection limit for <7nl ■ Following [50, [64J , 
we consider a (nearly spherical) survey of volume V. 
Assuming the Fourier modes are still uncorrelated and 
Gaussian distributed, the total signal-to-noise squared 
reads 



dominates the signal. When fc m j n <C fcmax, we can further 
simplify (S/N) 2 to 



i , io >+4^ Dr.) 



10- 



V 



/j- 3 Gpc 3 



4/3 



(34) 



Note the strong sensitivity of the signal-to-noise squared 
to the growth factor D(z) (For /nl non-Gaussianity, this 
dependence is only D(z)~ 2 ). 

To highlight the improvement one could achieve with 
future galaxy surveys, it is useful to first calculate the de- 
tection limit for the SDSS LRG sample centred at z ~ 0.3 
and covering a volume v ~ 2 /i~ 3 Gpc 3 . Assuming a lin- 
ear bias 6 = 2 and a skewness parameter ~2x 10~ 4 
appropriate for haloes of mass M ~ 10 12 — 10 13 M Q //i, 
the minimum ^nl detectable at the l-er level is ~ 10 6 
for a correction factor e K = 0.3 which we read off from 
Fig. [6j For a survey configuration analogous to SDSS- 
III/BOSS 2 , with central redshift z = 0.5 and a comov- 
ing volume V = 6 /i~ 3 Gpc 3 , the minimum <7nl would 
be ~ 4 x 10 5 for galaxies tracing haloes of similar mass 
and bias. Finally, for a configuration like EUCLID 3 
with a V = 100 /i~ 3 Gpc 3 survey centred at z = 1.4, 
the detection limit would be ~ 2.1 x 10 4 . Clearly, these 
limits are only indicative: they may be significantly im- 
proved by selecting highly biased, high redshift (single- 
or multi-)tracers. Nevertheless, this shows that future 
galaxy surveys should furnish interesting constraints on 
the size of the cubic coupling <?nl</> 3 - 



V 
4^2 



dkk 2 



Ab, 



(31) 



in the limit where sampling variance dominates the er- 
rors. Here, fc m ; n ~ tt/V 1 ^ is the smallest wavemode 
accessible and A: max is not necessarily finite since the in- 
tegral does converge as one takes fc max to infinity. Sub- 
stituting the expression Eq. (|23[) for the scale-dependent 
bias Ab K (k, <7nl) and setting T(k) = 1 over the wavenum- 
ber range across which the integral is performed, we ar- 
rive at 



1 



1 



(32) 



where 
k\ ~ 5.0 x 10-^g NL e 



2 - r -'> - 1:J - (1 }J h) [ ^— 1 //-Mpc 
D 2 (z) \ IO" 4 j 1 



(33) 

We have also assumed |<7nl| ^ 10 5 , such that \k 2 \ is at 
most of the order of k 2 lin and the term linear in Ab K /b 



C. Predictions for CMB temperature anisotropies 

The CMB trispectrum provides an alternative probe of 
local, non-quadratic correction to the Gaussian curvature 
perturbations, so it is interesting to assess the sensitivity 
of this statistics to the nonlinear parameter <?nl ■ 

The temperature anisotropy field is conveniently 
decomposed into spherical harmonics, AT(h)/T = 
J2hn a r^z m (fi)- A s shown in [2!j[3(|, statistical isotropy 
and invariance under parity transformation fi — * — n im- 
plies that the 4-point correlation of the spherical har- 
monic coefficients aY 1 takes the form 
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3 http: / /sci.esa.int/science-e/www/objcct/indcx.cfm?fobjectid=42266 
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Here, Q^(L) is the angular average trispectrum and 
brackets are Wigner-3j symbols. Statistical homogene- 
ity also implies that Q^^(L) is independent of position. 

The connected part of the trispectrum, Tf^(L), encodes 
information about non-Gaussianity and is obtained by 
subtracting a Gaussian piece constructed from the power 
spectra C;. Eq. (|3"5")) can be inverted with the aid of the 
orthogonality of the Wigner-3j symbols to form an esti- 
mator for the CMB trispectrum. 

The signal-to-noise for the CMB trispectrum T^j^L) 
summed up to a certain Z max is 30 
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(<C ^niax) 



£ £ 
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when cosmic variance dominates the errors. Otherwise, 
one shall include a contribution from the power spec- 
trum of the detector noise to the C; . Galactic foreground 
subtraction on a fraction 1 — / s k y would further reduce 
(S/N) 2 by a factor of / sky . 

Neglecting the ISW effect, the Sachs- Wolfe provides a 
useful order-of-magnitude estimate of the signal-to-noise 
as long as l max does not exceed < 100 0, H M [H[ . The 
calculation is performed in Appendix [Bl We find that the 
signal-to-noise can be recast into the compact form 



S Yr i ^ 9 2 ,2 f 1 r + \ 

]y ) (< 'max) = TJ^NlA, j g J ^ dx 



dxr?(x)s 2 (x) (37) 



where the auxiliary functions ri(x), si(x) and ti(x) are 
defined as 



n 



(x)=J2(2k + l)P k (x) 



k=2 
I 



i 

ti{x) = ^2(2k + l)k(k + l)P k (a 



(38) 
(39) 
(40) 



k=2 



Here, Pi{x) are Legendre polynomials. Note that we have 
excluded the monopole and dipole from the summation 
since these modes are unobservable. We have also as- 
sumed a nearly scale-invariant spectrum n s w 1. 

Figure. [5]shows the signal-to-noise ratio squared in the 
Sachs- Wolfe approximation for the CMB trispectrum as- 
suming 9nl = 1- Although this approximation breaks 
down for Z max > 100, we have extended the calculation 
up to Z max = 200 so as to extrapolate more robustly the 
^max-dependence to small angular resolution. A power- 
law fit to (S/N) 2 in the range 50 < Z max < 200 gives 
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FIG. 9: Signal-to-noise ratio squared for the CMB trispec- 
trum as a function of the maximum multipole Z max . We have 
assumed jnl = 1 and / sky = 1. 



Our results appear consistent with the findings of |30j 
shown in their Fig. 2. However, our constant of pro- 
portionality is about 20-30 times larger, presumably be- 
cause they adopted a lower fluctuation amplitude (com- 
pare also their prediction for the /nl model with that of 
[3ll]). Adding the information encoded in temperature- 
polarization trispectra may enhance (S/N) 2 by a factor 
of a few HI]. 

Assuming the scaling Eq. (|4"Tj) persists well beyond the 
range over which the Sachs- Wolfe effect dominates, the 
minimum <7nl detectable at l-er level is <?nl — 20, 7.9, 
3.2, 1.9 and 1.3 x 10 4 for Z max = 250, 500, 1000, 1500 
and 2000. A more realistic calculation should include 
the full radiation transfer function, detector noise etc. 
In this respect, detailed calculations have shown that, 
for the quadratic coupling /nl</> 2 , (S/N) 2 of the CMB 
bispectrum and trispectrum closely follows the behaviour 
obtained in the Sachs- Wolfe approximation [30|, |3j|. It 
seems reasonable, then, to expect that this is also true 
for 5nl^ 3 - 

While our predictions are qualitative they show that, 
for the WMAP CMB temperature measurement 4 (which 
we approximate as a noise-free experiment with ^ max ~ 
250), no detection of a significant trispectrum implies 
|.<?nl| < 2 x 10 5 at the l-er level. This is of the same order 
as the limit we derived from the QSO sample analyzed by 
For a PLANCK-like experiment 5 (Z max ~ 1500), no 
evidence for a trispectrum would imply |<?nl| < 1-3 x 10 4 
at the l-<7 level. This is comparable to the detection 



c) 



2.43xHT 17 ^ L 



10" 



2.6 



(41) 



4 http://map.gsfc.nasa.gov/ 

5 http:/ /sci.esa.int/science-e/www/area/index.cfm?fareaid=17 
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limit that could be achieved with an all-sky survey such 
as EUCLID. 



at the first order in the nonlinear parameters /nl and 
5nl- 



V. EFFECT OF NON-GAUSSIANITY WITH 
NON-ZERO / NL AND # NL 

In this Section, we examine the halo multiplicity func- 
tion and large-scale bias in numerical simulations of 
structure formation with non-zero coupling parameters 
(/nl, 9nl) = (±100, -3 x 10 5 ). We show that the results 
are consistent with those obtained from the simulations 
with non-vanishing <?nl solely. 



A. Mass function 

It is straightforward to calculate the fractional devia- 
tion from the Gaussian mass function, Eq. (|18|) . to non- 
zero /nl and <?nl- Again, we start with the MVJ formula 
and neglect second order corrections such as (crS 3 ) 2 etc. 
Adjusting the coefficient of the terms vcrSs and v 2 a 2 Si 
to that of the small v expansion obtained by [73| . we 
arrive at 



R(v, NL) = exp 



v v d(aS 3 ) 



(42) 



-—a Si 



v 2 d{a 2 S 4 ] 
4! d In v 



where the shorthand notation NL designates the combi- 
nation (/nl,5nl)- In Fig. [101 this theoretical prediction 
is compared R(v, NL) measured in non-Gaussian simu- 
lations of (/nl,5nl) = (±100,-3 x 10 5 ). We account 
for the fact that the amplitude of density fluctuations is 
renormalised by Sa§ ~ 0.0045. Fig. [TOl demonstrates that 
our approximation is in good agreement with the data, 
although it slightly overestimates the effect at z = 2 when 
/nl = —100. For /nl = 100, the positive and negative 
contributions from the quadratic and cubic coupling, re- 
spectively, almost cancel each other and flatten the frac- 
tional deviation over most of the mass range probed by 
the simulations. 

The scale-independent bias shift which arises from the 
change in the mean number density of haloes can again 
be estimated using the peak-background split. We find 



Abi(NL) 



1 



- 3v) a 2 Si + — {v 2 — 1) aS 



1 

+3' 

+2isS 



8v 



d{a 2 Si) vd 2 {a 2 Si) 



dlnu 
d{aS 3 ) 
rfln^ 



4! din v 2 
ld 2 (aS 3 ) 



6 d In v 2 



(43) 



B. Bias 

Having checked that the amount of non-Gaussianity 
in the mass function is also consistent with our simple 
theoretical expectation when both /nl and c/nl are non- 
zero, we now turn to the clustering of dark matter haloes. 
As shown in Appendix^ the non-Gaussian correction to 
the halo power spectrum can be written down as 



AP hh (fc) = ^bl5 c (z)a M (k)P^k) 



.9NL + ^/nL 



x ^(z)S^\M)a M (k)P^k) 



(44) 



If we set (?nl = and keep only the first two terms in the 
right-hand side, then the non-Gaussian (Eulerian) halo 
power spectrum can be cast into the form 

P hh (fc) = [b(M) + h h b^(k)] 2 P s (k) (45) 
where the scale-dependent bias parameter &^(fc) is 



fy(fc) =2[b(M)-l]6 c (z)a^(k) 



(46) 



Note that reference [64J obtained this relation by con- 
sidering the halo power spectrum implied by a bias re- 
lation that is a local mapping of the density field. In 
practice, the term proportional to P,p(k) is negligible 
as it contributes only at very small wavenumber k < 
0.001 /i _1 Mpc. The third term in the right-hand side 
of Eq. (|44p is derived in this paper for the first time. In 
the case c/nl = 0, its magnitude relative to the term 
linear in /nl is (l/5)f-^i J 8c(z)S 3 1 \M), which is approx- 
imately ~ 0.03 at redshift z = 1.8 and for a mass scale 
M = 10 13 M //i. Although its contribution becomes in- 
creasingly important at higher redshift, it is fairly small 
for the values of /nl considered here. Consequently, we 
shall employ the approximation 



A6(fc,NL) = e K Afe K (fc, 5NL ) + A& K (fc,/ NL ) 
+ [A6i(NL)+ei] . 



(47) 



to describe the non-Gaussian bias of dark matter haloes. 

The quadratic coupling /nl<^ 2 also affect the matter 
power spectrum at leading order [H, Hi[ , positive values 
of /nl increasing the small scale power. However, the 
relative size of this fc-dependent correction, /3 m (fc, /nl), 
is at a per cent level only in the weakly nonlinear regime 
k < 0.1 ftMpc -1 [53|,[85[ and fades rapidly as one goes to 
larger scales. We will thus neglect it in what follows. 

In Fig. [ITJ the result of measuring ratios of auto- and 
cross-power spectra in the simulations with (/nl,5nl) 
: ■ 1 00. ,'i x 10 5 ) is shown at < z < 1.5 for the haloes 
with M > 2 x 10 13 M Q /h. We do not quote error bars 
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v=6j k z)/o(U) 

FIG. 10: Fractional correction to the Gaussian multiplicity 
function of dark matter haloes as a function of the peak 
height v(M, z) for /nl = ±100 and a cubic coupling param- 
eter <;nl = —3 x 10 5 . The solid and dashed curves show the 
theoretical prediction Eq. (|42fl at z — and 2, respectively. 
Error bars denote Poisson errors. 



since the data points are obtained by averaging over two 
realisations only. The solid and dashed curves show the 
theoretical prediction Eq. (|4"?| . The value of the mul- 
tiplicative factor e K (b, <?nl) was obtained from the four- 
parameter formula Eq. (j28|) . while A&i + e\ was individ- 
ually fitted for each halo sample over the wavenumbcr 
range 0.005 < k < 0.03 /iMpc -1 . As can be seen, the the- 
oretical expectation Eq. (jl?) agrees reasonably well with 
the numerical data. This demonstrates that the range 
of validity of the non-Gaussian bias formula Eq. ([27]) ex- 
tends to smaller values of <7nl as well as models with 
non- vanishing /nl and <7nl- 

The lowest order, fc-dependent corrections to the Gaus- 
sian bias induced by the quadratic and the cubic coupling 
are fully degenerated in the halo power spectrum as they 
both scale as a^f(k) oc k^ 2 T(k)^ 1 . For the values of /nl 
and <7nl and the halo mass range considered here, the 
ratio A6 K (fc, <7NL)/A6 K (fc, /nl) increases approximately 
from 0.25 to 0.5 when the redshift increases from z — 
to 2. It is unclear whether higher order corrections could 
help breaking such a degeneracy. A more promising al- 
ternative may be to measure the bispectrum of dark mat- 
ter haloes, which carries much more information about 
the shape of the primordial three-point function than the 
power spectrum of bias tracers [40l. l4lj. However, this is 
beyond the scope of this paper. 



VI. DISCUSSION 

In this paper we explored the effect of a local cubic cou- 
pling 5nl^ 3 on the mass function and bias of dark matter 
haloes. We derived analytical expressions for these statis- 
tics and tested them against the outcome of numerical 
simulations. 

We showed that current theoretical predictions of the 
non-Gaussian correction to the mass function reasonably 
agree with the simulations. The LV formula [73| appears 
to provide a better fit to the data than the MVJ formula 
[74j , in agreement with some of the literature on the sub- 
ject [4(1 l8o| . The two approximations can be combined to 
provide an accurate description if one adjusts the low-v 
expansion of the latter so as to match that of the former. 

We found that the magnitude of the non-Gaussian 
scale-dependent bias Ab K (k, (?nl) is suppressed relative to 
a theoretical prediction based on the statistics of highly 
over dense regions, even on linear scales k < 0.01 /iMpc -1 . 
This suppression is stronger for the low biased sam- 
ples b < 3 and, at fixed value of 6, for positive values 
of <?nl- We were able to fit the measured halo bias 
at the expense of introducing two free parameters, e K 
and ei, that depend mostly on the halo bias b(M) and 
the coupling parameter c/nl- These parameters quantify 
the departure from the theoretical scale-dependent and 
scale-independent non-Gaussian bias correction, respec- 
tively. We provide a simple fitting formula for e K (b, <7nl), 
Eq.([2"8"]). which should be used when analyzing observa- 
tional data. In non-Gaussian simulations of the /nl type, 
the data also hint at a (possibly /NL-dependent) suppres- 
sion of the non-Gaussian scale-dependent bias relative to 
theory for wavemodes k < 0.03 /iMpc -1 [El, [II [El] , but 
the effect is much weaker than seen in our simulations 
of <7nl models. Clearly, these results require a better 
theoretical modelling of the non-Gaussian halo bias. 

Reference [EE1 argued that both the MVJ and LV an- 
alytic formula can be reconciled with measurements of 
the non-Gaussian fractional correction to the mass func- 
tion once non-spherical collapse is included. In practice, 
the critical density for collapse is replaced by <5 C — * ^/qS c , 
where the value q = 0.75 is obtained from a fit to the 
mass function measured in simulations (87(. Reference 
||| claimed that such a relation is a consequence of 
the diffusive nature of the critical threshold for collapse. 
Their model predicts q ~ 0.8, in good agreement with 
the findings of [El]. However, we found that substitut- 
ing 5 C — > i/qS c in Eq. (|18|) only modestly improve the 
agreement with the data. Regarding the non-Gaussian 
bias, it is not obvious how one could justify the replace- 
ment 5 C — > ^/qS c given that the linear bias of our (Gaus- 
sian) halo samples converges towards the spherical col- 
lapse prediction v 2 /S c for large peak height. 

A important ingredient is the choice of the halo iden- 
tification algorithm. While we used a spherical overden- 
sity (SO) finder, reference [H[ considered a Friends-of- 
Friends (FoF) finder with a linking length b = 0.2. The 
question of how the spherical overdensity masses can be 
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FIG. 11: Non-Gaussian bias correction measured in the simulations with (/nl,5nl) = (±100,-3 x 10 s ) for haloes of mass 
M > 2 x 10 13 M Q /h. Error bars are not shown as the data points are averaged over two realisations solely. 



mapped onto friends-of-friends masses remains a matter 
of debate [89[. Clealy however, since the peak height de- 
pends on the halo mass M through the variance cr(M), 
any systematic difference will be reflected in the value of 
v associated to a specific halo sample. This will in turn 
affects the size of the fractional deviation from the Gaus- 
sian mass function at some specified peak height. The 
sensitivity of the non-Gaussian mass function and bias 
to the halo finder will be presented elsewhere. 

The observational bound on /nl inferred by UK from 
the clustering of a high redshift sample of quasars can 
be straightforwardly translated into a limit on c/nl since 
Ab K (k, (?nl) also scales as k~ 2 T(k)~ 1 at low wavenumber. 
We have obtained 

- 3.5 x 10 5 < 5 NL < +8.2 x 10 5 (95%). (48) 

These are the first limits derived on g^L- While they 
are too weak to provide interesting constraints on infla- 
tionary scenarios such as the curvaton model, future all- 
sky redshift surveys should improve them by a factor of 
~ 100. Future CMB observations, including PLANCK, 



should also improve the limits derived here by an order 
of magnitude. With these improvements we expect that, 
in the next decade, realistic models of cubic type non- 
Gaussianity will be tested with real observations. 
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APPENDIX A: NON-GAUSSIAN BIAS IN THE 
HIGH PEAK LIMIT 

In this Appendix, we detail the derivation of the scale- 
dependent bias correction induced by the <?nl coupling 
to the two-point correlation of dark matter haloes. We 
follow (47| and approximate the latter by the two-point 
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correlation ^hh(r) of regions of the smoothed density field where the expression for the trispectrum of primordial 
5m with a peak height v ^> 1. curvature perturbation 



1. case (?nl 7^ only 



In this case, only the four-point correlations of the den- 
sity field contribute at the first order. The leading-order 
correction to the correlation of tresholded regions thus 
reads as 



i^ (4) (xi,xi,xi,x 2 ) + ^ (4) (xi,xi,x 2 ,x 2 

The four-point correlations £( 4 ) are Fourier transform of 
the trispectrum of the density field, Ta(ki, k 2 , k 3 , k.4), 
with conservation of momentum enforced. Linearly ex- 
trapolating the density field to present epoch, the latter 
can be expressed as 

T 5 (ki,k 2 ,k 3 ,k 4 ) = ^JJa M (fci)^ T$(ki, k 2 , k 3 , k 4 ) , 

(A2) 



T$(ki,k 2 ,k 3 ,k4) = 6 5N L[i 3 0(fci)P (fc 2 )P (fc 3 ) + (cyclic)] 
follows straightforwardly from the Fourier mode relation 



(2^ ( /7^^( k i)^(k2)^(k-k 1 -k 2 ) 

(A4) 

Density and curvature perturbations are related through 
the Poisson equation, whose scale-dependence is reflected 
in the function o^/(fc) = otn(k, z = 0) defined in Eq. ([6]). 
Combining these relations gives 



£ (4) (xi,xi,xi,x 2 ) 
£ (4) (xi,xi,x 2 ,x 2 ) 



65 



XL 



n 



65 



XL 



n 



d 3 h 



OiMikijP^ki) 

aM(h)P<p(ki) 



1 + 3 



1 + 2 



r 



P^fcm) P^k 123 ) 
Ptih) P (fc 3 ) 



aAf(fci23)e lkl23 ' r 

a M (k 123 )e^ r 



(A5) 
(A6) 



where wc have defined k 



ij—l 



k; for 



shorthand convenience. Since we will examine the effect 
of non-Gaussianity on Fourier space statistics only, we 

I 



take the Fourier transform of the four-point functions. 
After some simplification, we arrive at 



y i d 3 r^ 4 )(x 1 ,x 1 ,x 1 ,x 2 )e- <k - 1 



' d 3 k\ f d 3 k 2 

63xLaM(fc)P>(fc) / ,r,_^ a M(ki)P<p(ki) / 7^33 ^(^2)^(^2) 



a M (|k + ki 2 |) 

■d 3 h 



(2tt) 3 k ' vv +/ (2tt) 3 
P (jk+k 12 j) ' 



y"d 3 re (4) (x 1 ,x 1 ,x 2 ,x 2 )e^ k r = 6.g NL J^La M (ki)a M (\^+ ^iDP^P^ + k^) J 



(2tt) 3 

x a M (|k + k 2 |)P^(fc 2 ) 



d 3 k 2 
(2tt) 3 

' P„(|k + k 2 |) P„(|k + k 2 |) ' 

p (|k + k x |) p^k 2 ) 



(A7) 



(A8) 



For realistic values of the spectral index (n s ~ 1), the hand side of the above equalities formally diverges when- 
products aM(|k + ki|)P0 (|k + ki |) appearing in the right- ever k + ki =0 due to the ultraviolet divergence of 
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P^(fe). To cure this problem, one can set P$(k) = 
for sufficiently small wavenumbers or excise a thin shell 
centred at wavenumber fej from the integral. In the large- 



J 



scale limit k <C ki, the ratio P<^(|k + k,*|)/P^,(A;) vanishes 
whereas P^(|k + kf|)/P^(fej) tends towards unity. In this 
case, the above expressions reduce to 



/d 3 r^ 4 )(x l! x 1 ,x 1 ,x 2 )e 



a 3 5NL( r 4 4 1) (Af)aM(fc)P0(fc) 

f (Pki fcPk2 

+ 65NL«A/(fc) / 77, ^ OlM{kl)P^(kl) l-T—rzaM{k2)P(t>{k2)aM{kl2)P < t,(k\2), (A9) 



(27T) 



(27T) 



r 



and 



Jd 3 r^ (xi, xi, x 2 , x 2 )e" 



ik-r 



24 ffNL 



d 3 k! 
(2+j 3 



a 2 M (ki)P$(ki) ■ (AlO) 



Only the first term in the right-hand side of Eq. (|A9|) is 
not well behaved in the limit k — ► where it becomes 
proportional to fc™ s_2 . The second scales as k 2 , while 
the Fourier transform of £( 4 )(xi,xi,X2,X2) asymptotes 
to a constant. A similar decomposition also arises in 
/nl models. For this type of non-Gaussianity, the first 
order correction is furnished by the three-point function 
£( 3 '(xi, Xi, X2), whose Fourier transform can be split into 
the familiar term 2/nlc 2 ckm(&)P0(&), and a second piece 
given by 



1 



/ d ki 

hha M (k) j—— a l [ {k l )Pl{k l ). 



(All) 



In both quadratic and cubic local non-Gaussianity, the 
term proportional to ajvf(fc) can be neglected since, at 
the pivot point k = k$, its magnitude relative to the 
term involving a M {k)P^{k) is only 0(0.01) andO(10" 6 ), 
respectively. Moreover, it decreases as one goes to larger 
scales. By contrast, it is not so obvious how to handle the 
term (|ATol . In the non-Gaussian halo power spectrum, 
this term would appear multiplied by z^ 4 /(4<r 4 ), 



ft 



d 3 h 
J2nf 



aUk^ih) 



(A12) 



M 



10 13 M Q //i 



().:!7. r 1 



The approximation (second line) holds for 10 13 < M < 
10 14 M & /h. For 3nl = 10 6 and fer. ^ 3, this can be much 
larger than the typical shot-noise correction applied to 



the halo power spectra we measure in the simulations 
(see below). For (?nl = — 10 6 , this will certainly pro- 
duce a halo power spectrum which is negative at suffi- 
ciently low wavenumber. It is plausible that higher order 
counter-terms in the expansion Eq. (|2ip renormalises its 
value. However, such a calculation is beyond the scope 
of this paper, so the simplest choice is to ignore this 
term in the following of the analysis. Hence, we can 
approximate the non-Gaussian correction APhh in the 
limit of long wavelength k <C 1 by the Fourier transform 
of ^e^Hxi.x^xa^a)^ 4 , Eq. 



2. case $nl and /nl non-zero 

In addition to the first order trispectrum induced by 
<?nl</> 3 , Eq. (|A3|) . the quadratic coupling /nl^ 2 generates 
the bispectrum at leading order 

BsCki.ka.kg) = 2f NL [P(k 1 )P(k 2 ) + (cyclic)] , (A13) 

and an additional, albeit second order, contribution to 
the trispectrum (9fl, [M, [H , 

T$ (ki , k 2 , k 3 , k 4 ) = 4/ 2 L [P (fc 13 )P (fe 3 )P (t> (ki) 

+11 permutations] . (A14) 



The bispectrum (|A13p induces a three-point contribution 
(z/ 3 /er 3 )£( 3 '(xi,Xi,X2) to the pow er spectrum of biased 
tracers which is calculated in [47] . Upon Fourier trans- 
formation, it reads as 



APhh(fc) = 4/ NL 6^ c (z)a M (fc)P (fc) 



(A15) 



We follow the steps outlined above to calculate the con- 
tribution from the second order trispectrum (|A14[) . Af- 
ter some algebra, the Fourier transform of the four-point 
correlations of the density field can be written down as 
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,xi,xi,x 2 je 



= 8f^a M (k)P^(k) 

\ P*fo) 



d 3 ki 



(2tt) 3 
P (|k + k 12 |) 



d 3 k 2 



a-M (k2)P<t>(k2)aM (|k + k i2 1) 



(2tt) 3 

P^lk + kil) P^(|k + ki|) 



P (|k + k x |)P4|k + k 12 |) P (|k + k x |)P (|k + kia|) 



P^{ki)P^{k 2 ) 



P0(fc)P (fc 2 ) 



(A16) 



and 



d 3 r^{ 



x 1 ,x 1 ,x 2 ,x 2 Je 



- ik-r 



4 /nl , {2n)3 



d 3 ki f d 3 k 2 

a M (fci)aM(|k + ki|)P () ;,(fci) / — a M {k 2 )a M {\^- + ^-2\)P<t>{k 2 ) 



P (|k + ki 2 |) 

+P (fc) 



1 + 2 



(27T) 

P (|k + k!|) P^(|k + ki|)P^(|k + k 2 |) 



P0(fcl)P (fc 2 ) 



x + 2 P i 0k+k 1 Q p (|k + kiDp^ik + k 2 |) 



+2P (|k + k x |) 



P (fci) P0(fci)P (fe) 

P (fcl2) , P (fel 2 )" 



In the large-scale limit fc — > 0, these expressions asymptote to 



Jd 3 r£ {A) (x 1 ,xi,xi,x 2 )e 



-ikr 



|d 3 r^ 4 )(x 1 ,x 1 ,x 2 ,x 2 ) e — « 16/^ L ^P (fc) + 16/2 L |lAa| f ( fcl ) ai f (fe)P^(fc 2 )P (fc 12 ) 



x [P (fci)+P (fc 2 )] 



(A17) 



4/5l a 4 4 1) (M)a M (fc)P (fc) + 8/^a M (*) j p4 (*i 

7^aM(fc2)a M (A;i2)[P0(fei)P (fc 2 ) + P (fci)P (fc 12 ) + P|(fci 2 )] (A18) 
d 3 /ci -,,,,[ d 3 k 2 



(A19) 



Ignoring the second piece in the right-hand side of Eqs 
(|A18p and (j A19|) , the non-Gaussian correction to the halo 
power spectrum is then easily recast as Eq. (PHI) . 



APPENDIX B: SIGNAL-TO-NOISE FOR THE 
CMB TRISPECTRUM 



The formalism for the CMB trispectrum has been es- 
tablished in [2^, [3(| • The invariance of the 4-point har- 
monic function of the CMB temperature anisotropy field 
under the 4! permutations of the coefficients a™ i imposes 

constraints on the CMB trispectrum T^* (L) which can 



h+h 



phh 



(L') 



(Bl) 



be enforced by defining 

T}£(L)=Pti:(L) + (2L+l)Y, (-1) 

L' - 

h h L 
Is I4 V 

where curly brackets are Wigner-6j symbols, 

p}£( L ) = & L ) + c-i)^ 1 ' &) + (-i) L+h+u 

xt\±(L) + (-l) L +^+ l H\±(L)^ (B2) 

and the reduced trispectrum tj^(L) is symmetric under 
the exchange of its upper and lower indices and fully 
characterises the model. 

The expansion coefficients a] n are related to the pri- 
mordial curvature perturbation $(k) through 



a, 



4ir(-i) 1 



d 3 k 
(2^ 



$(k) PTi (fc)y ; m *(k) 



(B3) 
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where g Tl (k) is the radiation transfer function. The re- 
duced trispectrum can be calculated from this relation 
once the four-point function T$(ki, k2, IC3, IC4) is speci- 
fied. For a local cubic coupling 5nl</> 3 , 



t\l\l(L)= I drr*(3 l2 {r)p h (r)h HLl h, 



x [ mi (r)(3 h (r) + P h {r) m3 {r)] (B4) 



with 



Mr) 

w(0 



and 



47T 

We also use the notation 



*fc 2 P (% T! (fc)j z (fcr) 

T Jo 

2 f°° 

- dkk gNhg Tl (k)ji(kr) , 
n Ja 



h h L 




n =J(2h + l) 



x (21 j 



(B5) 
(B6) 

(B7) 

(B8) 



Note that most of the contribution to t\ ^(L) comes from 
a small volume centred at the comoving distance r* to the 
surface of last scattering. 

The Sachs- Wolfe approximation g Tl (k) w — j;(fcr i )/3 
valid at low multipoles I <C 1 00 p rovides a useful order- 
of- magnitude estimate @, H3, HH, [Hj]. In this limit, we 
can approximate Hi( r ) as — 5NL^n( r ~ r *)/(3 r *) since we 
assume c/nl independent of wavenumber. Hence, the re- 
duced trispectrum simplifies to 

t\±(L) « 9. 9NL Cf 2 w Cf 4 w (CT + Qt W ) ^^A^ , 

(B9) 

Inserting this expression successively into eqs (|B2[) and 
(IB1|) . the CMB trispectrum eventually reads as 
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T/ 1 / 2 (L) = — ffNL (2L + 1) 7T, 2 , , , 

3 4 27T 'l'a'sU 



(B10) 



/1 h L 


isw^sw^sw 



^3 h L 




+ (cyclic) 



where 



dkk 2 P^(k)j 2 (kr+) 



2irA 4> 
91(1 + 1) 



(Bll) 



The last equality assumes a nearly scale-invariant spec- 
trum n s ~ 1. The following relation between the Wigner- 
3j and 6j symbols (e.g., Appendix A of [29|). 



^(2^ + l)(-l) s+i3 -' 3 - mi - 



'1 ^2 '3 

W ^2 ^3 



'3 



h l[ 

m<i m'i — m 3 

h h h 
mi ™ 2 - m 



h 1' 3 
mi 7713 

h l'i 



(B12) 



where £ = ?i + £2 + 1[ + ^ ancl the value of m is set by the 
triangle condition, can be useful to derive Eq. (|B10[) . The 
signal-to-noise summed up to multipole Z max , Eq. (|36| . 
then becomes 



3\ 2 t , s_ f^V 2 

I (< Wj - I J 5NL ^ "iii a l,i« 

7 V 7 h>h>l 3 >U 



E n li 



[Cg w Cg w Cg w + (cyclic)] 1 

/l t2 £3 £4 



(B13) 



E + D 



L=0 



! 2 L\Vi 3 ? 4 L^ J 
I 



We can recast the sum over the diagonal modes L into 
a manifestly symmetric form with the aid of the Gaunt 
integral 



J dx P h (x)Pi 2 (x)P h (x) = ( q q q 
and the orthogonality relation 



(B14) 



^(2fc + l)P fe (x)P fc (2/) = 2S D (x - y) , (B15) 



k=0 



where 5d is the Dirac delta. We find 



E( 2i 



L=0 



Zl *2 i 





Z4 L 





(B16) 



dxP h (x)Pi 2 (x)P h (x)P u (x) 



There is a strict equality because the Wigner-3j symbols 
vanish for L > Eq. (|3"T|) for the signal-to- noise then 

follows by replacing the above equality into Eq. (|B13|1 and 
summing over all the 4! permutations of the quadruplet 
(hi hi h, h)- Although Eq. (|3T|) becomes computation- 
ally expensive when Z max 3> 100 (because we are sum- 
ming over redundant configurations), we found that it is 
quite efficient for Z max < 200. 

Following [3l|, we can roughly estimate the depen- 
dence of the signal-to-noise squared on Z max by con- 
sidering only the contribution of the L — 1 mode in 
Eq. (|B13[ ). Consequently, the product of the Wigner- 
3j symbols squared reduces to ~ Z1/3 5; 1 _i i ; 2 <5; 3 _i ! ; 4 and 
yields (S/N) 2 oc l^ax- However, including all L modes as 
in Eq. (|3"T|) gives a steeper dependence, (S/N) 2 cx l 2 ^^ 
(see Fig 12]), due to the fact that the Wigner-3j symbols 
decay slowly with increasing L. This is quite apparent in 
the classical limit l\, h, L ^> 1, where 



h h L 


1 



(B17) 



L[ {h + i 2 y„tfy ll2 [L 2 -(h-i 2 ) 2 ] 

7T J L J 



-1/2 
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Clearly, the terms in the summation Eq. (|B13[) decay only 
as 1/L for L ^> Z 1; Z 2 - By contrast, the second order 
contribution to the CMB trispectrum induced by the 
quadratic coupling /nl</> 2 adds an additional multiplica- 
tive factor of (Cf w ) 2 in the summation over the L modes 



which increases the relative contribution of the low-L 
modes (since these now decay as L~ 5 ). This is the reason 
why considering only L < 10 modes as done in [3l| still 
provides a good approximation to the signal-to-noise of 
the CMB trispectrum for the /nl model. 
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